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The dielectric and ferroelectric behaviors of a ferroelectric are substantially determined by its domain 
structure and domain wall dynamics at mesoscopic level. A relationship between the domain walls and high 
frequency mesoscopic dielectric response is highly appreciated for high frequency applications of 
ferroelectrics. In this work we investigate the low electric field driven motion of 90 "-domain walls and the 
frequency-domain spectrum of dielectric permittivity in normally strained ferroelectric lattice using the 
phase-field simulations. It is revealed that, the high-frequency dielectric permittivity is spatially 
inhomogeneous and reaches the highest value on the 90 "-domain walls. A tensile strain favors the parallel 
domains but suppresses the kinetics of the 90° domain wall motion driven by electric field, while the 
compressive strain results in the opposite behaviors. The physics underlying the wall motions and thus the 
dielectric response is associated with the long-range elastic energy. The major contribution to the dielectric 
response is from the polarization fluctuations on the 90° -domain walls, which are more mobile than those 
inside the domains. The relevance of the simulated results wth recent experiments is discussed. 

Substantial efforts have been devoted to ferroelectric (FE) materials and their applications in advanced 
electronic technologies 1,2 . The core ingredients in the physics of ferroelectrics include electric polarization 
P and its static/dynamic responses to electric field E ext or/and other stimuli. In the mesoscopic level, the 
dielectric and polarization responses of a ferroelectric are determined by the specific domain structure and 
realized by the domain wall motions. In particular, the high-frequency applications of ferroelectrics have received 
continuous attention. A number of FE materials in thin film form have been used for high-frequency devices in 
microwave and optical communications as well as computing applications because of their high dielectric 
permittivity 3,4 . A full understanding of the domain wall motion driven by either static or high-frequency dynamic 
force has been the core issue of the physics of ferroelectrics and the basis for their applications. 

The FE domain wall motion is a complicated process and depends on a series of intrinsic particulars such as 
defects, domain structures, and strains 5 7 . The consequence of strain, a topic to be dealt with in this work, has been 
highly concerned. Besides conventional scope of ferroelasticity, an externally imposed strain in a paraelectric 
surprisingly enables remarkable ferroelectricity 8 . Externally imposed strain can modulate the FE phase transitions 
too 9 " 11 . These strain effects raise particular attention to FE thin films deposited on rigid substrates, in which the 
induced strain can be controlled due to the lattice mismatch of the substrates with the films. It thus allows the 
performance improvement of the FE thin films by means of 'strain-engineering' the FE domains in the meso- 
scopic level. 

Along this line, the most interested issue is the strain effect in perovskite FE oxide thin films with tetragonal 
lattice symmetry like BaTi0 3 and PbTi0 3 . First, epitaxial BaTi0 3 and PbTi0 3 thin films deposited on substrates 
such as SrTi0 3 , LaA10 3 , and MgO, are often explored not only for fundamental understanding but also for 
potential applications, while textured polycrystalline films are also investigated. Second, these thin films offer 
regular twin-like (stripe-like) 90°-domain structure in coexistence with 180°-domains due to the intrinsic fer- 
roelastic effects. In most cases, the 90° ferroelastic domains are dominant, and the order parameter is strongly 
coupled to the strains. Furthermore, the substrate induced strain imposes a competitive or coherent coupling with 
the internal ferroelastic strain in the films, making the domain structure and domain wall motion complicated. 
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These behaviors, particularly in the high frequency range, are much 
less studied. It is believed that externally imposed strains have strong 
impacts on the domain wall motion and thus the high frequency 
dielectric permittivity 12-14 . 

A prominent feature of the strain effects is the 90° domain walls 
vibration in response to an ac electric field of high frequency, while 
the static responses have been well addressed. The dynamics of the 
domain wall may be characterized by the variation in dielectric per- 
mittivity as a function of the ac electric field (frequency w and ampli- 
tude £ 0 ), as long as the domain structure is well defined. A recent 
experiment demonstrating this dielectric response was carried out on 
clamped P^Zr^TiJOs (PZT) thin films deposited on Si sub- 
strates 15 . By creating cavities beneath the Pt/PZT/Si capacitors 
and cracking, these released PZT films show dramatic variation in 
the global dielectric nonlinearity and the frequency dependence as a 
function of mechanical clamping. The sequent piezoelectric 
examinations show that the increased local mobility of domain walls 
must be responsible for these behaviors, where "mobility" refers to 
the capability of domain walls in response to stimuli. This work 
among many others raises an important issue on how the 90°- 
domain structures in FE thin films couple dynamically with extern- 
ally imposed strains 16 . This question motivates us to investigate the 
dynamics of domain walls in the 90° domain structure modulated by 
strains, and consequently the dielectric responses at the mesoccopic 
level. 

While the strain coupling in polycrystalline PZT thin films may be 
too complicated, our purpose is to focus on tetragonal FE lattice (e.g. 
BaTi0 3 or PbTi0 3 ) and reveal the dynamics of the 90°-domain 
structure which couples the internal ferroelastic strains with extern- 
ally imposed normal strains, in response to electric field E ext . Instead 
of addressing the dynamic motion equation 1719 , we perform the 
phase-field simulations by which a continuous distribution function 
of the domain orientations is used to characterize the domain wall 
motion 20-21 . Such a phenomenological method enables us to directly 
examine the mesoscopic picture of these dynamic processes. We 
also investigate the underlying mechanism with which the strains 
modulate the domain wall motion in the energy landscape. In 
this phase-field model, the dipole-dipole and elastic interactions 
are considered 22 , since these long-range interactions play an import- 
ant role in configuring the domain structure 17-23 . The strain vector 
appears as an order parameter and can be modulated by external 
load. 

We consider a tetragonal FE system as approached by an L X L 
two-dimensional (2D) lattice with periodic boundary conditions 24 , 
and start from the Ginzburg-Landau theory. On each site, an electric 
dipole P(r) = (P x , P y ) normalized by a pre-factor P a 25 and an elastic 
displacement vector u(r) = (u x , u y ) are imposed. We clarify that our 
simulation is not necessarily associated with a thin film, although the 
electric dipoles in our simulation can only relax on the in-plane 
configuration. An extension of the results to other geometry is plaus- 
ible. For simplification, our simulation will not deal with the finite 
boundary problem since it will introduce complicated corrections to 
the mechanical balance 25 , while open-circuit boundary issue will only 
be discussed briefly. Given this assumption, the lattice can be 
mapped into a small region of real material. By these approximations, 
we can focus on the details of the local dielectric response in meso- 
scopic scale. 

The electric field E ext =(E x , E y ) is confined on the lattice plane. The 
total free energy for this FE lattice can be written as 26 : 



F = Fid + F g + Fdd- 



Pes H - *si 



(i) 



where F\d, F g , Fdd, F e i, F es , and F se are the Landau potential, gradient 
energy, dipole-dipole interaction, elastic energy, electrostrictive 
energy, and electrostatic energy, respectively. Term F id extending 
to the sixth-order is written as: 



Fu = A, [P 2 x +P 2 y ) +A n [Pt + Py) +AnP 2 x P 2 r 

<<•> ^ (2) 

+A m (p« + P 7 6 ) +An 2 (P 2 x P y +P x Pf 

where A\, A n , A 12 , A lu and A 112 are the Landau expansion coeffi- 
cients and A 1 =A 10 (r-r o ). The lowest-order expression of term F g i&: 



- G 44 {Px,y + Py,,) 2 + G 44 (P x ,y - P y , x ) 2 



(3) 



where Py=3P/3»"j and rj={x, y), and Gn, Gi 2 , G 44 and G' 44 are the 
gradient energy coefficients 20 . Terms Fdd and F se can be written 
respectively as: 



Fdd~- 



4ns 



EE 

<f> <j> 



p(r,yp(rj) 3P(n)-(n - rj) [PtyjHr, - rj)} 



-.(4) 



= (E x ,Ey), 



(5) 



where Fd d includes two contributions. One is the anisotropic part 
Fdep(Po) and the other is the isotropic part F'dd(8P) if we treat 
P(r)=P 0 +5P(r) where P 0 is the spatially homogeneous average 
polarization and r is the spatial vector. It is noted that P 0 can be 
defined by: 



SP(r)dV = 



dP(r)dr 3 = 0, 



(6) 



where V is the volume of lattice. 

Term Fdep(Po) describes the depolarization energy and shares the 
same form as F se by difference of a 1/2 factor. In real system, the 
depolarization field is compensated by free charges for the cases with 
open boundaries. For the present case, no free charge is included. The 
effective electric field is independent of the depolarization field and 
simply equivalent to E ext . In experiments, the electric field is usually 
applied via the electrodes in close-circuit which introduces free 
charge onto the interface. This will cancel out the depolarization 
field. Term F d d counts an integration over the whole lattice, and a 
realistic calculation is done either by the Fourier transformation or 
by finite truncation treatment 27-28 . For a 2D lattice, the finite trun- 
cation is a sufficiently accurate approximation as long as the trun- 
cating distance R is big (P=8 in our simulation) 24-29-30 . 

The elastic energy F el yields: 



Fd= E 1 Cn ( u Xx + u j,y) + C uu x , x u hy + - C44U*. 



(7) 



Ujj = duj I drj, Uij = 8ui/ 8rj + duj/dri if i 



where Cy is the elastic stiffness tensor which has only three inde- 
pendent elastic constants for a square lattice. 
The electrostrictive energy F es is: 

F es = - ^ u^iqnPl + quPp + u^/qnPy+quP^ + u^yq^Py, 

(8) 

where i?ii = C 11 Q 11 +2C 12 Qi2, <2i2=C 11 Q 1 2+2C 12 (Qii + Qi2), and 
1^44 — 2C44 Q 44 are the effective electrostrictive coefficients, Q« are 
the electrostrictive coefficients in external stress free state. 
Mathematically, we can also separate the total strain ;/y into a homo- 
genous component if 0 , and a heterogeneous one <5)7y which denotes 
the microscopic strain distribution at site (;', j) 22-25 : 



nii(r)=n 0 + 5nii{r)- 



(9) 
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here external strain r\ 0 is imposed either by a substrate or an external 
load. Similar to generally accepted treatment 20-22,28 , order parameter 
u{r) is reduced to a function of P(r) 31 . The step-by-step procedure of 
this derivation can be found in the Supplementary materials. To 
access the dielectric response over frequency domain, we compute 
dielectric permittivity e(cu) over a broad range of frequency. The 
algorithms are presented in the Method section. The parameters 
for the simulations are listed in Table I, and for the details of the 
simulation procedure one can refer to earlier works 32,33 . 

Results 

Domain wall kinetics under dc electric field. We first investigate the 
90°-domain wall motion driven by a static {dc) E ext . To characterize 
the wall motion, we define the domain width / deviating from the 
equilibrium width Z 0 under zero electric field. Thus, Al=l-l 0 stands 
for the offset in response to E ext . In order to minimize the effect of 
domain wall irregularity, parameters I and l 0 are calculated indirectly. 
A set of rectangle-like regions, whose two sides are on the domain 
walls to enable their areas as big as possible, are taken. We sum the 
areas of these rectangles. The local domain width I is obtained by 
dividing the area by its length and then performing statistics over 
sufficient number of rectangles in the same region. 

In Figure 1 is shown the kinetics of field-driven domain wall 
motion, given different strains from slightly compressive ones to 
tensile ones. Both E ext and t] 0 are along the y-axis, thus the I indicates 
the width of a 2 -domain where the dipoles align along they- axis. The 
evaluated Al(t) data are plotted in Figure 1(a) where the applied 
strains are labeled numerically. For each strain, Al(t) shows similar 
behavior and increases rapidly in the early stage and tends to be 
saturated at Al max in the late stage. Nevertheless, Al max depends on 
ri 0 , and the bigger r\ 0 the smaller Al max . In Figure 1(b) we plot the 
evaluated l o 0lo) f° r the a 2 -domain and corresponding Al max (ri 0 ), 
showing the monotonous increasing of Z 0 ('?o) an d decreasing of 
Al max (ii 0 ). For definition of each symbol and motion pattern of 
domain wall, one can refer to the Supplement materials. 

The kinetics of wall motion was investigated earlier 34 and can be 
described by a simple model 17,23 , yielding the kinetic equation 
Al(t) = a{l-e~ bt ) where a=Al max and b are the fitting parameters. 
As shown in Figure 1(c), the monotonous decrease and increase of 
a and b as a function of i] 0 respectively imply that the a 2 -domain 
extension becomes tougher upon the transition from compressive 
strain to tensile one. It is also possible to evaluate the initial motion 
speed of the a 2 -domain walls by taking v 0 =dAlldt\ t _, 0 =ab, as shown 
in Figure 1(d), consistent with the above argument. In other words, 
for compressive strain (;?o<0)> the motion speed of domain walls is 
higher than that for tensile strain (;/ 0 >0). 

For the cases with strain 1] 0 not parallel to E ext , the above analysis is 
qualitatively correct. The tensile strain makes the domains with P//;y 0 
wider and the compressive one makes it narrower. This behavior 
allows a competition between the strain effect and electric field effect, 
and they can be coherent or cancelled depending on their orientation 
relationship. However, when the strain contains shear components, 
the situation becomes more complicated and extensive calculation is 
not discussed here. 
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Figure 1 | Evaluated parameter \1 as a function of time t given a series of 
labeled numerically (a), i 0 and A/ max as a function of i/ 0 respectively (b), 
parameters a and 6 as a function of ij 0 respectively (c), and initial domain 
wall motion speed v 0 as a function of //„ (d). E ext =6\A i \P a . 

Domain wall vibrations and dielectric response: r|o = 0. Based on 
the above result, one understands that the 90°-domain structure 
clamped by external stress has different stability characteristics 
from the stress-free state. This difference can be discussed in the 
clamped 90°-domain structure driven by the ac-electric field, 
characterized by variation of dielectric permittivity as a function of 
r\ 0 . In particular, the dielectric response in the frequency domain is 
related to the wall motion. 

We first address the dielectric response in ?/ 0 = 0- Figure 2(a) and 
(b) show the real part s r {f) and imaginary part £,(/) at three 0 angles. 
Here a small E 0 =Q.6\A 1 \P a is chosen, with the field direction defined 
by angle 0 between the x-axis and E ext . In general, £ r (f) decreases 
gradually with increasing / for all the three cases, while £,(/) shows 
two peaks at characteristic frequencies^ ~ 0.1t~' andf H ~ 7t _1 , 
which are respectively referred as the low-/ and high-/ dispersions. 
Here t= l/jA^D is the characteristic time for electric dipole flip (see 
the Method section for details). The striking feature is the low-/ 
dispersion anisotropy, i.e. the ^-dependence which is the most 
remarkable at 0 = 90° (and 0=0 too, four-fold symmetry). This in- 
dependence is shown by e r (0) at f=0.01z' 1 as an example in 
Figure 2(c), characterized by the typical periodic variation. 

The above behaviors can be understood by investigating the 
instant evolution of domain structure. By turning E ext from 0=0 to 
0= 180°, one checks the domain evolution and dielectric dispersion. 
The low-/ dispersion is related to the domain wall vibrations, while 
the high-/ dispersion is attributed to the flip of individual dipoles. In 
fact, the dispersion peak around f H is also observed in the mono- 
domain lattice. 

Now we investigate the origin for the low-/ dispersion anisotropy 
in the domain scale. In f;o = 0> the spatial distributions of s r at 
/=0.05t _1 in four 0 angles are presented in Figure 3(a)~(d), respect- 
ively, where the color scales the intensity (e r =0.0~1.0). It is imme- 
diately seen that the dielectric permittivity mainly comes from the 
contribution of wall vibrations, while those dipoles deeply inside the 



Table 1 Physical parameters chosen for the sirru 


lation (t '= A] D) [24, 25]. All these parameters appear in the 


dimensionless form 


Parameter (unit) 


Value 


Parameter (unit) 


Value 


Parameter (unit) 


Value 


1 

A? 2 (A^Po/lA,!) 
Gf, (G n /A 2 |A, |) 
G'44 (GWAtKI) 

C44 (C44/|A,|Po 2 ) 

giLilW |Ai |) 


64-256 
2.50 
1.60 
0.80 
0.543 
0.0157 


A* ( A, |) 
A, n (AniPo/lA, 
GT 2 (G 12 /A 2 |A, 
C,, [C u /\A,\Po) 
9,, (qii/|A, |) 

T* (t) 


-1.00 
) 0.49 
0.00 
2.75 
0.143 
0.0004 


At, (A u Py\A,\) 
An 2 (A112P0/IA1 1) 
Gl 4 (G 44 /A 1 2 |A 1 |) 
C* 2 [C }2 /\A,\Po] 
q 12 (qi 2/I Ail) 


-0.24 

1.20 

0.80 

1.79 
-0.0074 
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Figure 2 | Simulated dielectric permittivity spectra: (a) real part and (b) 
imaginary part over a wide range of frequency, at 0=90°, 1 12°, and 135°, 
respectively, (c) Dielectric permittivity real part ?, r as a function of angle 9 
at a constant /=0.01t~'. 

domains contribute little. This characteristic makes the dielectric 
response very specific. 

At 0=45°, the domain structure and dielectric response are shown 
in Figure 3(a). Again, the dielectric response mainly comes from the 
contribution of electric dipoles on the walls and near regions. This 
feature can be understood from the electrostatic energy -P-E ext . The 
walls have the highest local mobility. No matter E ext is positive or 
negative, there is always one of the neighboring two domains, inside 
which all the dipoles take the 135° angle from E ext . This domain will 
shrink while the other will expand, making the wall move easily. This 
is also the reason why £ r is the highest at 0=45°, as seen in Figure 2(c). 

At 0=90° (or 0=0), as shown in Figure 3(b), the wall mobility is 
slightly lower, and e, r is lower too. In this case, E ext is parallel to the 
dipoles in one domain and perpendicular to those in the other. At 
0=112° and 0=135°, as shown in Figure 3(c) and (d), respectively, 
the wall mobility falls down even further, resulting in even lower E r . 
The wall mobility becomes the lowest at 0 = 135°. Figure 3(d) shows 



y/% 




0.5 




Figure 3 | Snapshoted patterns of dielectric permittivity real part s r at 
constant frequency /= 0.05 r -1 with angle #=45° (a), 90° (b), 112° (c),and 
135° (d). 



that almost the whole lattice has the identical £ r except the very thin 
and dim lines on the walls. The reason is that term -P-E ext in two 
neighboring domains are almost equivalent and they compete with 
each other, hindering the wall motion. The dielectric response is 
weak with no frequency dispersion. 

We finally check the dielectric response in the high-/ range and 
one example is given in Figure 4(a) and (b) at/=6t _1 for two specific 
0 angles: 0 = 90° and 0=135°. The dielectric distribution over the 
whole lattice is roughly homogeneous. For details, one looks at the 
case in Figure 4(a) and finds only weak color contrast between 
the two neighboring domains a 1 and a 2 . The domain a l whose 
dipoles align along the x-axis but perpendicular to E ext shows slightly 
higher e r than that of domain a 2 . The reason is obvious that the 
dipoles in domain a 1 are more fluctuating than those in domain a 2 . 

Domain wall vibrations and dielectric response: ti o >0. Now we 

discuss the cases with iio^O. We only present in details the results on 
normally strained lattice. When the lattice is normally strained, the 
domain structure is deformed. To clarify the similarity and difference 
between the strain-free and strained lattices, we consider the simplest 
situation: )? 0 =0.7% under E ext with /=0.05t~ 1 and £0=0.6^!?,,, 
both aligned along the y-axis (0 = 90°). We present in Figure 4(c) 
and (d) the strain-free domain structure and strained structure 
respectively, as well as the spatial distribution of e r . For the two 
cases, the distributions are similar in amplitude but the high-e r 
spatial profiles across the domain walls are much wider for the 
strain-free lattice than those for the strained lattice. In the quali- 
tative sense, the dielectric permittivity averaged over the whole 
strained lattice is lower than that over the strain-free lattice, at 
least in the low-/ range. This (j 0 -dependence can be further 
illustrated in Figure 5(a), where £ r (f) and £,(/) at three tensile 
strains (>/ 0 = 0> 0.2%, 0.7%) are plotted. Both the real and imaginary 
parts are remarkably suppressed by tensile strain, and the f L has a 
slight shift towards the high-/ side. We also calculate the instant 
responses of the Al and polarization component P y against E ext at 
the three tensile strains, shown in Figure 5(b). The one-to-one 
correspondence between Al, P y , and E cxt , respectively, is observed. 
The responses of Al and P y are synchronous with E ext in the low-/ 
range, and delayed in the high-/ range. Furthermore, the wall 
vibration amplitude is suppressed by the tensile strain. 
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Figure 4 | Snapshoted patterns of dielectric permittivity real part e r at 
constant frequency/=6r _1 with angle 0=90° (a) and 135° (b). At constant 
frequency /=0.05t~' with different tensile strains along the y-axis: (c) ij o = 0 
and (d) ); 0 = 0.7%. The red lines below the snapshots indicate the alignment 
of s r along the [1,1]. Hereafter, £0=0.6^1?,,. 
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Figure 5 | (a) Simulated dielectric permittivity spectrum around 
/=0.05t~' given i?o=0, 0.2%, and 0.7% along the y-axis. (b) From top to 
bottom: parameter A/as a function of time t given );o=0, 0.2%, and 0.7% at 
/=0.05t _1 , the y component of total polarization as a function of time f, 
and the ac electric field E ext as a function of time t. The ac electric filed is 
along y-axis, E 0 =0.6\A 1 \P a . 

The above results refer to the simplest situation. When r\ 0 and E ext 
are aligned along arbitrary directions independently, more complex- 
ity is seen in terms of the domain structure evolution and dielectric 
response. For some specific geometry, the strain and electric field 
may compete with each other. Some more discussion will be given 
below. However, in general, the calculated results are qualitatively 
similar to those for the simplest situation: the tensile strain sup- 
presses the domain wall vibration, thus reducing the dielectric 
permittivity. 

Domain wall vibrations and dielectric response: n 0 <0. Now we 

check the cases with fjo<0. Referring to the results under static (dc) 
E ext , as shown in Figure 1, one sees that the static compressive strain 
assists the 90° domain wall motion. It is thus expected that the 
dielectric permittivity in compressed lattice will increase. The r\ 0 - 
dependences of £ r (f) in three compressive strains (;/ 0 =0, —0.1%, 
— 0.15%) along the y-axis are plotted in Figure 6(a), consistent 
with the expected results. Similar behaviors are observed for the 
strain and electric field both applied along the x-axis. 

We also check the situations where E ext is not parallel to 1] 0 . One 
example is shown in Figure 6(b) for ;; 0 >0 and Figure 6(c) for f; o <0, 
where E ext is along the y-axis and i] 0 is along the x-axis. Our extensive 
calculations establish the qualitatively similar behaviors: the tensile 
strain suppresses the 90° domain wall motion and thus the dielectric 
permittivity, while the compressive strain enhances the domain wall 
motion and the dielectric permittivity, no matter whether E ext is not 
parallel to i-] 0 or not. We summarize spectrum e r (/, >? 0 ) in Figure 7(a) 
and (b). At/</ L , such as/=0.0lT _1 and 0.05t~\ e r falls gradually with 
increasing rj 0 from r] 0 <0 to f] 0 >Q, while this tendency becomes neg- 
ligible as f>fi such as/=0.5T~' since the single dipole response 
becomes dominant at this frequency. The overall evolution of e r (/, 
i]o) is plotted in Figure 7(b). 




Figure 6 | (a) Simulated dielectric permittivity spectrum around 
/=0.05t _1 given );o=0, —0.1%, and —0.15% along the y-axis with the ac 
electric along the y-axis too. The calculated real part e r at tensile strain (b) 
and compressive strain (c) perpendicular to the ac electric field along the y- 
axis. £ 0 =0.6lAilP fl . 

Discussion 

It should be mentioned that the periodic boundary conditions allows 
the total strains in all directions to be exactly balanced out, and the 
strain effect is of long-range and coupled with mechanical boundary 
conditions. However, if other mechanical boundary conditions are 
considered, such as free boundary, the total strains can be partially 
relaxed through the free boundaries. The 90°-domain walls can be 
more moveable and thus show more significant response to electric 
field, as partially discussed in Ref. 15. 

To understand the simulated results, one may give additional 
discussion by looking at the energy landscape and comparing the 
simulated results and experiments, although relevant experimental 
data are really rare. 

Energy landscape. We calculate the elastic energy distribution 
associated with the 90° domain structures. For a single domain 
lattice, the total strain field is homogeneous, thus accommodating 
extremely large elastic energy. The single domain is decomposed into 
the 90° domain structure so that the total elastic energy can be 
relaxed. Owing to the lattice volume conservation, for one domain, 
the compressive strain along one direction (e.g. the x-axis) is always 
accompanied with the tensile strain along the other direction (the y- 
axis), and verse vice. One can refer to the Supplement materials for 
the strain distribution. 

First, we consider the t] 0 =0 case. Given a static electric field along 
the y-axis, the domain walls move into the a Y -domain in compensa- 
tion with the extension of the a 2 -domain width. The continuous 
shrinking of the (^-domain is accompanied with the increasing mag- 
nitude of elastic strain (e xx and e yy ) inside the (^-domain. Since the 
total elastic energy is proportional to if, the rapidly enhanced total 
elastic energy inside the shrunk ^-domain acts as the resistant force 
against the electric field, responsible for eventual termination of the 
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Figure 7 | (a) Simulated dielectric permittivity real part £ r as a function of 
jj 0 with the ac electric field frequency /=0.01t _1 , 0.05t _1 , and 0.5t _1 . (b) A 
2D plot of £ r as a function of tj 0 andf. The ac electric field and strains are all 
along the y-axis. E 0 =0.6\A 1 \P a . 
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wall motion when the ^-domain becomes sufficiently narrow, as 
shown in the left column of Figure 8(a) and Figure 8(b), where the 
total elastic energy V e \{x, y) is plotted. This explains why Al tends to 
be saturated at Al max as time goes infinitive. 

This scenario applies to the lattice with externally imposed strain. 
Take ;; 0 =0.4% along they-axis as an example. At f=0 with E ext =0, 
the whole elastic energy distribution on the right column of Figure 8 
shifts upward with respect to the case of 170= 0. However, one can 
observe a much higher energy distribution in the shrinking a x - 
domains. Therefore, the t^-domains are already under highly tensile 
state along the x-axis even at E ext =0, when an externally imposed 
strain along the y-axis is applied. In this case, a static electric field 
along the y-axis may further force the wall moving into the d\- 
domain but will be highly resisted by the high strain energy, as shown 
in the right column of Figure 8. Therefore, additional shrinking of the 
fl r domain driven by the electric field would become even tougher, 
given that r] o =0A% already makes the ^-domain narrow. 

If strain r\ 0 is applied along other orientations, no matter whether it 
is parallel to E ext or not, the wall motion behaviors remain qualita- 
tively similar. One example is given by applying E ext along the x-axis 
where E ext and r\ 0 are normal to each other. As shown in Figure 8(c), a 
strain r\ 0 along they-axis enhances remarkably the total elastic energy 
in the a 2 -domain, making the further shrinking of the a 2 -domain 
more difficult than the case with ?7 0 = 0. Therefore, the above-dis- 
cussed results don't lose the generality. 

Comparison with experiments. The implication of the above 
simulation data can be checked by comparing the simulated results 
with experimental data on the contribution of domain wall 
vibrations to the dielectric response in BaTi0 3 (BTO) and 
Pb(Zr 1 . x Ti x )0 3 (PZT) ceramics 35 37 , although these data have not 
yet received confirm from the direct detecting of the mesoscopic 
scale dielectric distribution across the 90° domain walls. Therefore, 
such comparisons may not be so direct in quantitative sense. For 
PZT, experiments 35 showed that above 10 11 Hz the real part of 
dielectric constant begins to drop, due to the frozen dipole 
oscillators themselves. In our simulation, the corresponding 
frequency is assigned as f H = 7t~ 1 in Figure 2(a). Thus we have 



Figure 8 | Evaluated spatial contours of elastic energy F eI in the 90°- 
domained lattice. The left column (a) refers to ij o =0 and the right one (b) 
refers to >7 0 =0.4%. The top row refers to E ext =0, the middle row refers to 
t = 720t after E ext ( dc) along the y-axis applies to the lattice, and the bottom 
row refers to f=720r given E M (dc) along the x-axis. E ext =6\A 1 \P a . 



-10" 



and the 



the characteristic inverse time scale t 
calculated f L ~10' J Hz. 

Regarding the PZT ceramics, experimental data 36 also indicate the 
dielectric dispersion around 1.0 GHz, which was argued to arise 
from the vibrations of the frozen 90° domain walls. For the dielectric 
permittivity magnitude, we take the dimensionless factor 
(|Ai|£*)-'~10 3 , where A x is 3.8(T-479) X 10 5 CT 2 m 2 N taken from 
PZT and £* is the vacuum permittivity. The calculated real part of the 
dielectric permittivity at 0=45° is —200, as shown in Figure 2(c), and 
the maximum of the imaginary part is —100. Indeed, experimental 
measurements 37 on ceramics PZT (x= 0.48-0.52) gave the difference 
in the real part between the extremely high frequency and the very 
low frequency, which is 250-500. The peak of the imaginary part is 
150-265. These values agree roughly with our simulation results 
here. The above comparisons allow us to claim that the present 
calculations are reliable even in quantitative sense. 

We also find some experimental results about the effects of tensile 
and compressive strains on the overall dielectric constant. In micro- 
wave frequency range, it was shown that dielectric constant and 
tunability of BTO films grown on MgO gradually decrease as the 
in-plane strain goes from the compressive type to tensile type 38 . 
Similar behaviors were observed in polycrystalline Ba 0 .6Sr 0 .4TiO 3 
thin films upon increasing tensile strain 39 . Those results are consist- 
ent with our calculation in a qualitative sense. Unfortunately, none of 
those experiments establishes a clear logic between the dielectric 
response and the 90° domain wall motion, although the frequency 
range is properly related to the 90° domain wall vibration. 
Establishing this logic is a technical challenge for experimentalists. 
To reveal the domain wall response to external electric field, one may 
map the local polarization response in real-space. Nevertheless, this 
task becomes difficult for such a high frequency. In this sense, the 
present simulation seems to be unique for bridging the relationship 
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between the dielectric response in microwave frequency and the 
strain via the microscopic domain scale. 

Methods 

Calculation of ac dielectric permittivity. The temporal evolution of the dipole lattice 
is tracked by solving the Ginzburg-Landau (TDGL) equation which takes the 
following form: 

dP(r,t) _ SF 



a 



SP(r,t) ' 



(10) 



where t is time scaled in unit x with z ~ l = \Ai\D, and D is the kinetic coefficient. We 
also calculate the dielectric permittivity as a function of E ext . For a dielectric system 
that cannot polarize instantaneously in response to an electric field, the total electric 
polarization as a function of t can be described as: 



p( t )= 



B(t-t')E exl (t')dt\ 



(11) 



where P(t) is a convolution of electric field E(t) at previous times with time-dependent 
permittivity e(f)=£ r (f) + i£i(f)' Therefore, its Fourier transformation can be directly 
written as 

P(oj) = s(co)E exl (co), (12) 

where P(co) and E ext (co) are the Fourier transformations of P(t) and E ext (t) 
respectively 32 . For the case of ac sine electric field such as E ext (t)—E 0 -sin{co Q t) with 
CQ 0 = 2nf 0 and co=2nf E cxt (oj) can be directly written as: 

E ext (co)=EoS(co — coo) 3 (13) 

where E' 0 is the coefficient of Fourier transformation and co 0 is the frequency. 
Polarization P{co) can be calculated by Fourier- transforming the temporal evolution 
spectrum P(r, t) in Eq.(10). In details, the real-spatial spectrum of P(r, t) at site r and 
its spatial average P(t) over sufficient number of time periods is calculated by solving 
Eq.(10). The Fourier-transformed frequency spectrum P(r, co) can be expressed as: 



P(r,0)) = \ 
Jo 



= P{r,t)e- ia>t dt, 



(14) 



which is used to transform P(r, t) and P(t) into frequency domain to obtain P(r, co) 
and P(oj). Subsequently, one can compute the corresponding dielectric permittivity 
e(r, co 0 ) at site r and its spatial average e(co 0 ) at frequency co Q using Eq.(12). 
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